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Abstract: We compute the expansion of the 3-d Lattice QCD free energy to four-loop 
order by means of Numerical Stochastic Perturbation Theory. The first and second order 
are already known and are correctly reproduced. The third- and fourth-order coefficients 
are new results. The known logarithmic divergence in the fourth order is correctly iden- 
tified. We comment on the relevance of our computation in the context of dimensionally 
^ ■ reduced finite temperature QCD. 
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1. Introduction 



The QCD free energy density (or the pressure of the quark-gluon plasma) is a good observ- 
able to study the deconfinement phase transition [Qj. The goal is to study the transition 
between the realm of low-temperature hadronic matter, where confinement is the main 
physical phenomenon, and the quark-gluon plasma phase that is realized at high tempera- 
tures, which in turn is governed by asymptotic freedom. In the latter phase the pressure is 
given by the Stefan-Boltzmann limit of an ideal gas of non-interacting particles, p oc T . 
Ideally, one would like to undertake lattice simulations across the phase transition up 
temperatures at which the pressure exhibits a purely perturbative behavior. In practice, 
however, the convergence properties of the perturbative expansion are poor at temperatures 
which are not asymptotically large Q, while on the other hand computational resources 
limit the highest temperatures at which lattice simulations can be performed (a fair limit 
is some 4-^5 times the transition temperature T c ~ 200MeV). 

Dimensional reduction || is a strategy to fill the gap one is facing, and in fact has 
been applied to the problem in question The setup is as follows. One starts with the 
full theory (Ad QCD) and as a first step matches this to a 3dSU(3) gauge theory coupled 
to a Higgs field in the adjoint representation. This theory can then be matched to 3d pure 
gauge SU(3), which captures the ultrasoft degrees of freedom. Both these reductions have 
been successfully performed in a continuum (perturbative) scheme, i.e. MS. 3d pure gauge 
SU (3) then has to be treated non-perturbatively, the only practical method being lattice 
measurements. 

In order to incorporate these lattice measurements into the reduction setup, it is essen- 
tial to know the relation between the two regularization schemes. This is the point where 
Lattice Perturbation Theory (LPT) comes into play. Due to the superrenormalizable na- 
ture of the 3d theory, all divergences can be computed perturbatively. This allows a clean 
matching of the schemes in the continuum. Computing at high orders in LPT is not a 
simple task (in the present case we need g s order - note that this means four loops for the 
free energy, but three loops for the plaquette), and that is why we make use of Numerical 
Stochastic Perturbation Theory (NSPT) §]. 

We recall the definition of the free energy density / 

Z = jDUe~ Sw[u] =e"^ / , (1.1) 

where the (Wilson) pure gauge action reads 

Sw = Po - lip) , (1.2) 
p 

with (3q = 2N C / '(a 4_d <7o) denoting the the standard dimensionless (bare) lattice inverse 
coupling in d dimensions, while lip is the basic plaquette 

lip = -^(Tr Up) , (1.3) 
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which is to be computed at any point on any independent plane according to 



Up = U^{n) = U^n)U u {n + M )t/t( n + v)U}(n) . (1.4) 
To compute the free energy one can now revert to the computation of the plaquette 

(l-Hp) = Z- 1 j DUe~ Sw M (I -Up) (1.5) 



2 a d d , „ 2 a d d A ,, 
lnZ = "773 TTnS-ta/) • 



d(d-l)V5/9b 1) dfi T 

Hence, given a weak-coupling expansion of the plaquette 

(1 - nP> = ^^ + ^- + ^- + ^- + - (L6) 
it follows that 

2a d ,1 ^ (AT , Ar ... „ c 2 (iV c ,d) c 3 (iV c ,rf) c 4 (iV c ,aQ 

^- iy (-/) = c (iV c ,rf) + c 1 (iV c ,d)ln/? ^ ^ ... . (1.7) 

We now specialize to QCD (iV c = 3) in d = 3 dimensions, where g§ ~ a -1 and hence 
A) = 6/(a^o)- The previous formula reads (from here on, q = Ci(AT c = 3, d = 3)) 

| i/ = a~ 3 (co + Cl ln/3 ) - a" 2 | ffo 2 - a" 1 ^ - ^ 5o 6 + 0(a) . (1.8) 

In order to control the matching to continuum one then needs the first four coefficients 
in the expansion of the basic plaquette. Note however that it is already known from a 
computation in the continuum that at four loop level there is a logarithmic infrared (IR) 
divergence ||. One of the aims of our computation is to recover the scheme-independent 
coefficient of this logarithm, while fixing the lattice constant which is left over once an IR 
regulator has been chosen. In Eq. ( |1.8|) we have put a tilde on C4 to denote that the IR 
divergence has to be isolated and subtracted in a convenient scheme. Later on, the lattice 
size L will act as the IR regulator. Going back to the lattice coefficients themselves, the 
first and the second ones are already known |7j . We will give them in Section ^. The third 
and the fourth ones are the goal of the present work, a task which one can manage within 
our computational scheme 1 . 



2. Computational setup 

Our computational tool is NSPT. Computing to /3q~ 4 order requires to expand the field up 
to /3 -4 order ||, that is 

8 

C/» = l + ^/3 ^C7«(n). (2.1) 
i=i 

1 In 4d the expansion of the plaquette is known via NSPT up to a much higher order. It is interesting 
to compare the two different situations. In 3d there is the additional subtlety of the IR divergence. On 
the other hand, the dimensionful nature of the coupling in 3d makes it possible to single out the different 
divergent contributions in Perturbation Theory. The situation is much more involved in 4d (see ||). 
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We write the expansion in terms of the Un field. One could also express everything in terms 

of the Lie algebra field A^(n) = Ya=i A) ^ i n )-> the relation being U^n) = exp(A M (n)). 
Whichever one uses, one should keep in mind that perturbation theory amounts in any case 
to decompactify the formulation of lattice gauge theory. The expansion in terms of the 
ujfi (n) is easier to manage from the point of view of computer data organization. Eq. ( |2.1| ) 
is the expansion to be inserted in the Langevin equation 

d t U v = [-iVS[U v ] - i V ]U v , (2.2) 

r\ being a Gaussian noise. The equation has to be integrated in a convenient (time) dis- 
cretization scheme. Our choice is the Euler scheme as it was proposed in |J. This amounts 
to introducing a time step e. As usual, the solution is recovered by working in a region 
where the time step corrections are linear (Euler scheme is a first order scheme) and extra- 
polating to e — > 0. We computed the expansion in Eq. ( |1.6| ) for (1 — lip) on different lattice 
sizes ranging from L = 5 to L = 16 (up to three loops we also performed the computations 
on a L = 18 lattice) 



, c {L) c {L) c (L) c {L) 

(l - n ^WW-- <2 - 3) 



We then extrapolated the infinite lattice size results according to 



(L) . di 

1 1 + L3 



ln fc L 



<r = ^+££ d r ] -^f- (<=2,3) (2.4) 

j=ji k=0 



=cf) In L 3 + C4 + jr^dP 

3=34, k=0 



ln k L 
Li 



These asymptotic forms are basically dictated by Symanzik's analysis [|H]]. In particular 
they include the contribution from subleading logarithms (they are suppressed by inverse 
powers of L). For each i the index k runs up to i — 1, which equals the number of loops 
(remember that in terms of the plaquette we are computing up to three loops, i.e. i = 1 
is the tree level). The index j counts the subleading contributions coming from inverse 
powers of L. As it appears from the last line of the previous formula, the finite volume 
also acts as the IR regulator needed at order (this is instead a leading logarithm). 
Some comments are in order at this point. The final errors on infinite volume results are 
dominated by this extrapolation process. 

Trying to assess the effect of the subleading logarithms, it turns out that both the 
range of our data and our statistical errors do not allow to distinguish between a logarithm 
and a constant. Hence we will only give (effective) extrapolations based on pure power- like 
fits. The spread of the results comes from the indetermination on the (inverse) powers to be 
included in the fit. This is not surprising, since in NSPT there is no control on what in the 
language of Feynman diagrams would be contributions coming from different diagrams {i.e. 
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(£>) 
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(L) Jln), T 3 
c\ — c\ In L 


5 


2.5455(13) 


1.8682(45) 


CT nnn / rjp\ 

5.990(26) 


o r nn do) 

25.99(18) 


oi 00/1 o\ 

21.28(18) 


6 


2.6536(8) 


1 C> (\( ' 0/01 \ 

1.8968(31) 


6.200(19) 


27.66(14) 


22.41(14) 


7 


2.6580(8) 


1.9095(30) 


6.307(21) 


28.68(15) 


oo no/i p* \ 

22.98(15) 


8 


2.6615(6) 


1.9226(23) 


6.408(16) 


29.66(14) 


oo p" T / 1 a\ 

23.57(14) 


9 


2.6630(6) 


i nooo ( o o \ 

1.9288(22) 


£i A O A f 1 0\ 

6.484(18) 


on a a { ~\ n\ 

30.44(16) 


o A nn / 'X n\ 

24.00(16) 


i n 

10 


O C f O O / A\ 

2.6638(4) 


1.9340(15) 


6.519(13) 


on ni / 1 o\ 

30.91(13) 


O/l 1 p/i o\ 

24.16(13) 


11 


O G£2 AtL ( A\ 

2.6645(4) 


1.9381(14) 


6.574(11) 


31.53(14) 


O/l C 1 / 1 A \ 

24.51(14) 


12 


2.6650(3) 


1.9413(12) 


6.591(11) 


31.67(15) 


24.39(15) 


13 


2.6653(3) 


1.9423(12) 


6.621(11) 


32.27(18) 


24.76(18) 


14 


2 6656(3) 


1 9436(12) 


6 288(11) 


32 37(16) 


24 64(16) 


15 


2.6662(2) 


1.9455(10) 


6.652(10) 


32.84(19) 


24.91(19) 


16 


2.6657(2) 


1.9442(8) 


6.658(9) 


33.28(19) 


25.16(19) 


18 


2.6663(2) 


1.9489(7) 


6.715(8) 






oo 


2.6666(1) 


1.955(2) 






25.8(4) 



Table 1: The coefficients c\ at the various lattice sizes and their infinite volume extrapolations 
Ci. For the last order we report both and — InL 3 ; the latter is the quantity to be 
extrapolated. For the error on c 3 see text and Figure |2[ 

sums). We then try different choices of the powers and then compare the corresponding 
X 2 's. This process does not select a definite set of powers: the better choices (see Figure |2|) 
turn out to be comparable with respect to the resulting x 2 - The quoted values for the c, 
(i > 1) together with the associated errors embrace the range of the outcomes. 

Most of the computer simulations have been performed on a PC cluster the Parma 
group installed one year ago. This is made of ten bi-processor Athlon MP2200. A pro- 
gramming environment for NSPT for Lattice Gauge Theory was set up in C++. This was 
in part inspired by the TAO codes we use on the APE machines and for a large fraction 
based on the use of (C++ specific) classes and methods to handle lattice and algebraic 
structures. Needless to say, this part of the work will be useful in other applications of 
NSPT. The results we report come out of 6 months of runs on the above mentioned cluster. 
Some more statistics came from another PC cluster more recently installed in Parma. The 
latter is a blade system based on 14 Intel Xeon 2.0 GHz. 

3. Results 

In Table [l] we present the results we obtained for the coefficients at various values of 
L. As in Eq. ( |l,6j ) and Eq. ( J2.3| ), we only give the coefficients of order /3^ n : the coefficients 
odd in go (i.e. of order (3 ^ n+1 )/ 2 ) we re verified to be zero within errors. In the last line 
one can read the values extrapolated to infinite volume. For the fourth order we present 
both the bare coefficients cjj^ and the subtracted ones — InL 3 ; the latter is the 
finite quantity one is interested in at L = oo. 
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In Figure |] we plot the values of the coefficients at the various lattice sizes together 
with the interpolating finite size corrections. Again, for the fourth order we plot both the 
bare coefficients and the subtracted ones — In L 3 . An obvious benchmark for 
our computations is the first order, whose value is = 8/3 * (1 — 1/L 3 ). In this (trivial) 
case one knows the result both at finite and at infinite volume. That is why in Figure [| 
we also plot the known finite size corrections for c±. 

Another benchmark is the second coefficient, which is also found in agreement with 
the diagrammatic studies in J7J. From Ref. [|ll]] one reads C2 = 1.9486, however without an 
error estimate. 

As it was already pointed out, for orders higher than the trivial one, the quoted 
errors of the infinite volume-extrapolated values are dominated by the form of the fitting 
polynomials in Eq. (2.4). Still, the final errors are acceptable. Note the asymmetric error 
for C3. The more conservative lower bound takes into account a choice for the subleading 
powers of L which results in a worse x 2 > see Figure |2[ It is interesting to compare the 
result C3 (our first original result) with the one conjectured in jjll]] from the hypothesis of 
the dominance of a given contribution: C3 ~ 7.02. This conjecture turned out to be not 
too far from the result. 
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Figure 2: Different fits for C3 . The upper 
curves takes 1/L as the leading power, while 
the lower (dashed) one takes 1/L 2 . Note that 
the latter results in a higher \ 2 (0-7 vs 0.55). 
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Figure 3: Different fits for c. 
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In these fits the value for c\ is fixed to 
be the analytically known one. Within the 
range of the data the fits almost coincide. 



Let us now discuss the IR divergence at order {3q 4 - As is well known, it is difficult to 
fit a logarithm. Still, we obtain enough evidence for it. By this we mean the following. One 
can take different approaches to the fit of the last line of Eq. Q2.4| ). One possibility is to 
include no logarithmic correction at all. A second one is to include a logarithmic correction 
whose prefactor is a fitting parameter. A third possibility is to include a logarithmic 
correction whose prefactor cf n ^ equals the result which has already been obtained in the 
continuum computation of ||: cjf = 81 (688 — 1577r 2 /4)/(47r) 4 = 0.9765. By varying 
the choice of the inverse powers included in the fit, the first case (no log) yields values 
of x 2 which are systematically worse than in the other two cases: as expected, the data 
appear to prefer the inclusion of the logarithm. By fitting both C4 and we obtained 
C4 = 24.5(2.0) and = 1.1(2), a result fully consistent with ||. In the third case (see 
Figure |3|) we obtain C4 = 25.8(4), getting a smaller error like expected. For the meaning 
of the quoted errors, see the discussion in Section |2[ 

4. Conclusions and perspectives 

We computed the first four coefficients in the expansion Eq. ( |1.8|) of the plaquette in 
3d pure gauge SU(3) theory, from which one can trivially obtain the expansion of the 
free energy at four loops. For the first two coefficients the already known results have 
been correctly reproduced. The third coefficient (the first original result of this paper) is 
connected with the mildest power divergence to be subtracted from simulations data for a 
lattice determination of the free energy. 

The fourth coefficient was known to be logarithmically divergent, a result that was 
reproduced. Let us further comment on this point. A priori one can use any IR regulator 
in order to extract a finite part for the four loop contribution: the finite volume (the one 
we used in the present work), a mass (a very popular IR cutoff in the continuum) or the 
coupling itself (since it is dimensionful in 3d) . Notice that the latter is in a sense the natural 
choice for computer simulations. Obviously each choice defines a scheme of its own. While 
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the coefficient of the logarithm is universal, there are of course specific constants relating 
the different schemes. As already pointed out, the coupling itself is the most natural 
regulator for computer simulations, even if there is no simple way to take it as the cutoff 
in perturbation theory. Ultimately, we are interested in the matching between the lattice 
and a continuum perturbative scheme (to be definite, MS). The idea is to take the same 
IR regulator both in continuum and in lattice perturbation theory, which most naturally 
would be a common mass for all tree-level propagators. Since the same mismatch will be 
present in both computations with respect to the data coming from computer simulations, 
that mismatch will cancel in the matching. Employing massive propagators in NSPT will 
therefore be the natural extension of the approach presented here. 
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